﻿Imports DotSpatial.Topology, DotSpatial.Data, DotSpatial.Symbology, System.Threading.Tasks
#Const SingleThreaded = False

''' <summary>
''' Encapsulates routines to convert polygon or point shapefile to raster, determine IDW weights, and reassign grid values based on attribute values
''' </summary>
''' <remarks></remarks>
Public Class clsGrid
    Public Enum enumCellSizeOption
        CellSize
        NumCells
    End Enum

    Private lyrSource As FeatureLayer = Nothing
    Private lyrClip As FeatureLayer = Nothing
    Private FieldName As String = ""
    Private FileName As String = ""
    Private CellSizeOption As enumCellSizeOption = enumCellSizeOption.CellSize
    Private CellSize As Double = 50
    Private NumClosestPoints As Integer = 5
    Private MaxDistance As Double = 500
    Private NumCats As Integer = 50
    Private IDWExponent As Integer = 1
    Private lstNeighbors(0) As List(Of Tuple(Of Integer, Double)) 'indices of nearest points in feature layer and their weight (1 to NumClosestPoints)
    Friend lyrGrid As FeatureLayer = Nothing

    ''' <summary>
    ''' Instantiate class
    ''' </summary>
    ''' <param name="lyrSource">Layer with points or polygons having numeric values in attribute field</param>
    ''' <param name="FieldName">Name of attribute field in source</param>
    ''' <param name="FileName">File name to use for new grid file (polygon shape file)</param>
    ''' <param name="CellSizeOption">How is CellSize value to be interpreted</param>
    ''' <param name="CellSize">Size of rectangular cell (map units) or number or columns or rows to create</param>
    ''' <param name="NumClosestPoints">Number of closest points to use for interpolation</param>
    ''' <param name="MaxDistance">Points farther away than this distance are ignored</param>
    ''' <param name="IDWExponent">Inverse distance weighting exponent</param>
    ''' <param name="NumCats">Number of categories to create for symbology scheme</param>
    ''' <param name="lyrClip">Layer to use for clipping of grid</param>
    Sub New(lyrSource As FeatureLayer, FieldName As String, FileName As String, CellSizeOption As enumCellSizeOption, CellSize As Double, NumClosestPoints As Integer, MaxDistance As Double, IDWExponent As Integer, NumCats As Integer, lyrClip As FeatureLayer)
        Me.lyrSource = lyrSource
        Me.FieldName = FieldName
        Me.FileName = FileName
        Me.CellSizeOption = CellSizeOption
        Me.CellSize = CellSize
        Me.NumClosestPoints = NumClosestPoints
        Me.MaxDistance = MaxDistance
        Me.IDWExponent = IDWExponent
        Me.NumCats = NumCats
        Me.lyrClip = lyrClip
    End Sub

    Public Function CreateGrid() As FeatureLayer
        Try

            Dim fsGrid As New FeatureSet(FeatureType.Polygon)
            fsGrid.DataTable.Columns.Add("Result", GetType(Double))

            Dim starttime As DateTime = Now

            With lyrClip.Extent

                If .Height = 0 Then Return Nothing

                'determine number of rows and columns, accounting for aspect ratio of clipping extents

                Dim NumRows, NumColumns As Integer
                Dim dxdy As Double = .Width / .Height
                If dxdy > 1 Then 'is wider than tall
                    If CellSizeOption = enumCellSizeOption.NumCells Then
                        NumColumns = CellSize + 1
                        NumRows = NumColumns / dxdy
                    Else
                        NumColumns = Math.Truncate(.Width / CellSize) + 1
                        NumRows = NumColumns / dxdy
                    End If
                Else
                    If CellSizeOption = enumCellSizeOption.NumCells Then
                        NumRows = CellSize + 1
                        NumColumns = NumRows * dxdy
                    Else
                        NumRows = Math.Truncate(.Height / CellSize) + 1
                        NumColumns = NumRows * dxdy
                    End If
                End If

                Dim dx = .Width / (NumColumns - 1)
                Dim dy = .Height / (NumRows - 1)

                'create list of simplifies external boundary polygons so testing during clipping is faster
                Dim lstPoly As New List(Of Polygon)

                If TypeOf (lyrClip) Is PolygonLayer Then
                    For Each f In lyrClip.DataSet.Features
                        Dim poly = DouglasPeuckerReduction(CType(f.BasicGeometry, Polygon), Math.Min(dx, dy))
                        lstPoly.Add(poly)
                        'fsGrid.AddFeature(poly) 'temporary
                    Next
                End If

                For r = 0 To NumRows - 1
                    If Not GISForm.AnimationInfo.IsAnimating Then WriteProgress(GISForm, r, NumRows - 1, "Creating grid...")
                    Dim o As New Object
#If SingleThreaded Then
                    For c = 0 To NumColumns - 1
#Else
                    Parallel.For(0, NumColumns,
                                    Sub(c)
#End If
                                        Dim x = .MinX + c * dx
                                        Dim y = .MinY + r * dy

                                        Dim pt As New Point(x + dx / 2, y + dy / 2)
                                        Dim found As Boolean = False
                                        If TypeOf lyrClip Is PolygonLayer Then
                                            For Each poly In lstPoly
                                                If poly.Contains(pt) Then found = True : Exit For
                                            Next
                                        Else 'clipping to point layer, just create grid over entire extents and don't do clipping
                                            found = True
                                        End If

                                        If found Then
                                            Dim Coords As New List(Of Coordinate)
                                            Dim rand As New Random()
                                            For k = 1 To 4
                                                Coords.Add(New Coordinate(x + Choose(k, 0, 1, 1, 0, 0) * dx, y + Choose(k, 0, 0, 1, 1, 0) * dy + rand.Next(0.01)))
                                            Next
                                            SyncLock (o)
                                                With fsGrid.AddFeature(New Feature(FeatureType.Polygon, Coords))
                                                    .DataRow("Result") = DBNull.Value
                                                End With
                                            End SyncLock
                                        End If
#If SingleThreaded Then
                    Next c
#Else
                                    End Sub)
#End If
                Next r
            End With

            With fsGrid
                .Projection = New DotSpatial.Projections.ProjectionInfo
                .Projection = lyrSource.Projection
                .InitializeVertices()
                .UpdateExtent()
                .InvalidateVertices()
                .InvalidateEnvelope()
                .SaveAs(FileName, True)
            End With

            lyrGrid = GISForm.GetLayer(FileName)
            If lyrGrid IsNot Nothing Then GISForm.GISMap.Layers.Remove(lyrGrid)
            lyrGrid = GISForm.GISMap.AddLayer(FileName)

            If Not GISForm.AnimationInfo.IsAnimating Then WriteProgress(GISForm)

            Debug.Print("CreateGrid: " & Now.Subtract(starttime).TotalSeconds)

            Return Me.lyrGrid

        Catch ex As Exception

            ErrorMsg(, ex)
            Return Nothing

        End Try
    End Function

    Public Sub ComputeWeights()

        Try

            'todo: Consider using non-Euclidian interpolation (like ArcMap) where points that cross a boundary (e.g., clip layer) are not considered
            'todo: Weights are computed in advance; if doing animation and none of the closest 5 points exist, will assign NULL to grid cell

            If lyrGrid Is Nothing OrElse lyrSource.DataSet.NumRows = 0 Then Exit Sub

            ReDim lstNeighbors(lyrGrid.DataSet.Features.Count - 1)

            Dim MaxDist As Double = 0.0 'set max dist based on largest polygon

            Dim lstPoly As New List(Of Coordinate)
            For f = 0 To lyrSource.DataSet.Features.Count - 1
                With lyrSource.DataSet.Features(f)
                    If TypeOf lyrSource Is DotSpatial.Controls.MapPointLayer Then
                        lstPoly.Add(.Coordinates(0))
                    Else
                        lstPoly.Add(.Centroid.Coordinates(0))
                        MaxDist = Math.Max(MaxDistance, Math.Max(.Envelope.Height, .Envelope.Width))
                    End If
                End With
            Next

            If MaxDist = 0.0 Then MaxDist = Me.MaxDistance

            Dim starttime As DateTime = Now

#If SingleThreaded Then
        For g = 0 To lyrGrid.DataSet.Features.Count - 1
            'WriteProgress(GISForm, g, lyrGrid.DataSet.Features.Count - 1, "Computing weights...")
#Else
            Parallel.For(0, lyrGrid.DataSet.Features.Count, Sub(g As Integer)
#End If
                                                                Dim dictDist As New SortedDictionary(Of Double, Integer)
                                                                Dim pt = lyrGrid.DataSet.Features(g).Centroid.Coordinates(0)
                                                                Dim rand As New Random()
                                                                For f = 0 To lstPoly.Count - 1
                                                                    Dim delta As Double = (rand.NextDouble - 0.5) * 0.00001 'to avoid duplicate keys
                                                                    If dictDist.Count < NumClosestPoints Then
                                                                        dictDist.Add(pt.Distance(lstPoly(f)) + delta, f)
                                                                    Else
                                                                        If Math.Abs(lstPoly(f).X - pt.X) <= MaxDist AndAlso Math.Abs(lstPoly(f).Y - pt.Y) <= MaxDist Then
                                                                            Dim dist As Double = pt.Distance(lstPoly(f))
                                                                            If dist < MaxDist Then dictDist.Add(pt.Distance(lstPoly(f)) + delta, f)
                                                                        End If
                                                                    End If
                                                                Next

                                                                lstNeighbors(g) = New List(Of Tuple(Of Integer, Double))

                                                                With lstNeighbors(g)

                                                                    Dim sumwt As Double = 0.0

                                                                    For Each kv In dictDist
                                                                        If .Count < NumClosestPoints AndAlso kv.Key < MaxDist Then
                                                                            Dim wt As Double = 1.0 / kv.Key ^ Me.IDWExponent  'inverse distance raised to a power
                                                                            sumwt += wt
                                                                            .Add(Tuple.Create(kv.Value, wt)) 'add index and weight if closest one or within radius of influence
                                                                            If .Count = NumClosestPoints Then Exit For
                                                                        End If
                                                                    Next

                                                                    'recalculate the weights
                                                                    If sumwt > 0 Then
                                                                        For i = 0 To .Count - 1
                                                                            .Item(i) = Tuple.Create(.Item(i).Item1, .Item(i).Item2 / sumwt)
                                                                        Next
                                                                    End If
                                                                End With
#If SingleThreaded Then
        Next
#Else
                                                            End Sub)
#End If

#If SingleThreaded Then
        For g = 0 To lyrGrid.DataSet.Features.Count - 1
            'WriteProgress(GISForm, g, lyrGrid.DataSet.Features.Count - 1, "Computing weights...")
#Else
            Parallel.For(0, lyrGrid.DataSet.Features.Count, Sub(g As Integer)
#End If
                                                                Dim feat = lyrGrid.DataSet.Features(g)
                                                                Dim pt As New Coordinate(feat.Centroid.Coordinates(0).X, feat.Centroid.Coordinates(0).Y)
                                                                Dim dictDist As New SortedDictionary(Of Double, Integer)
                                                                Dim rand As New Random()
                                                                For f = 0 To lstPoly.Count - 1
                                                                    Dim delta As Double = (rand.NextDouble - 0.5) * 0.00001 'to avoid duplicate keys
                                                                    If dictDist.Count < NumClosestPoints Then
                                                                        dictDist.Add(pt.Distance(lstPoly(f)) + delta, f)
                                                                    Else
                                                                        If Math.Abs(lstPoly(f).X - pt.X) <= MaxDist AndAlso Math.Abs(lstPoly(f).Y - pt.Y) <= MaxDist Then
                                                                            Dim dist As Double = pt.Distance(lstPoly(f))
                                                                            If dist < MaxDist Then dictDist.Add(pt.Distance(lstPoly(f)) + delta, f)
                                                                        End If
                                                                    End If
                                                                Next

                                                                lstNeighbors(g) = New List(Of Tuple(Of Integer, Double))

                                                                With lstNeighbors(g)

                                                                    Dim sumwt As Double = 0.0

                                                                    For Each kv In dictDist
                                                                        If .Count < NumClosestPoints AndAlso kv.Key < MaxDist Then
                                                                            Dim wt As Double = 1.0 / kv.Key ^ Me.IDWExponent  'inverse distance raised to a power
                                                                            sumwt += wt
                                                                            .Add(Tuple.Create(kv.Value, wt)) 'add index and weight if closest one or within radius of influence
                                                                            If .Count = NumClosestPoints Then Exit For
                                                                        End If
                                                                    Next

                                                                    'recalculate the weights
                                                                    If sumwt > 0 Then
                                                                        For i = 0 To .Count - 1
                                                                            .Item(i) = Tuple.Create(.Item(i).Item1, .Item(i).Item2 / sumwt)
                                                                        Next
                                                                    End If
                                                                End With

#If SingleThreaded Then
        Next g
#Else
                                                            End Sub)
#End If
            'WriteProgress(GISForm)

            Debug.Print("ComputeWeights: " & Now.Subtract(starttime).TotalSeconds)


        Catch ex As Exception

            ErrorMsg(, ex)

        End Try

    End Sub

    Public Sub UpdateGrid(Optional FieldName As String = "")

        Try

            If lyrGrid Is Nothing OrElse lyrSource.DataSet.NumRows = 0 Then Exit Sub

            If FieldName <> "" Then Me.FieldName = FieldName Else FieldName = Me.FieldName

            With lyrGrid.DataSet.Features

                Dim starttime As DateTime = Now

                Dim o As New Object

#If SingleThreaded Then
            For g = 0 To .Count - 1
#Else
                Parallel.For(0, .Count, Sub(g)
#End If
                                            Dim val As Double = 0 'product of value and weight
                                            Dim sum As Double = 0 'sum of weighting factors

                                            If lstNeighbors(g) Is Nothing Then
                                                sum = Double.NaN
                                            Else
                                                For Each t In lstNeighbors(g) 'item 1 is index, 2 is weight
                                                    Dim f As Feature = lyrSource.DataSet.Features(t.Item1)
                                                    If IsNumeric(f.DataRow(FieldName)) Then
                                                        val += f.DataRow.Item(FieldName) * t.Item2
                                                        sum += t.Item2 'note: for some time steps, predetermined nearby stations may have no data; so only consider those that do
                                                    End If
                                                Next
                                            End If

                                            SyncLock (o)
                                                If sum = 0 Then 'no valid nearby points found
                                                    lyrGrid.DataSet.DataTable.Rows(g)("Result") = DBNull.Value
                                                Else
                                                    lyrGrid.DataSet.DataTable.Rows(g)("Result") = Math.Round(val / sum, 5) 'did not fill datatable until i used this approach????? 
                                                End If
                                            End SyncLock

#If SingleThreaded Then
            Next g
#Else
                                        End Sub)
#End If

                lyrGrid.DataSet.Save()
                Debug.Print("UpdateGrid: " & Now.Subtract(starttime).TotalSeconds)
            End With

        Catch ex As Exception

            ErrorMsg(, ex)

        End Try

    End Sub

    Friend Shared Sub ClipWithPolygon(lyrPoly As PolygonLayer, lyrClip As PolygonLayer)
        Try

            'create list of simplifies external boundary polygons so testing during clipping is faster
            Dim lstPoly As New List(Of Polygon)

            'get cell size from first polygon
            Dim dx As Double = lyrPoly.DataSet.Features(0).Envelope.Width
            Dim dy As Double = lyrPoly.DataSet.Features(0).Envelope.Height
            For Each f In lyrClip.DataSet.Features
                Dim poly = DouglasPeuckerReduction(CType(f.BasicGeometry, Polygon), Math.Min(dx, dy))
                lstPoly.Add(poly)
            Next

            With lyrPoly.DataSet
                For f = .Features.Count - 1 To 0 Step -1
                    Dim found = False
                    For Each poly In lstPoly
                        If poly.Contains(CType(.Features(f).BasicGeometry, Polygon)) Then
                            found = True
                            Exit For
                        End If
                    Next
                    If Not found Then 'delete the feature
                        .Features.RemoveAt(f)
                    End If
                Next
                .InitializeVertices()
                .UpdateExtent()
                .InvalidateVertices()
                .Save()
            End With


        Catch ex As Exception

            ErrorMsg(, ex)

        End Try

    End Sub

    ''' <summary>
    ''' Create coloring scheme for raster layer based on current values
    ''' </summary>
    Friend Sub CreateScheme(Optional MinValue As Double = Double.MaxValue, Optional MaxValue As Double = Double.MinValue)
        Try

            If lyrGrid Is Nothing Then Exit Sub

            With GISForm.AnimationInfo

                Dim starttime As DateTime = Now

                Dim ps As New PolygonScheme
                ps.Categories.Clear()
                ps.SuspendChangeEvent()
                ps.SuspendEvents()

                If MinValue = Double.MaxValue AndAlso MaxValue = Double.MinValue Then 'desired range not passed; determine from data in source layer

                    For Each feat In lyrGrid.DataSet.Features
                        Dim v = TestNull(feat.DataRow("Result"), 0)
                        MinValue = Math.Min(MinValue, v)
                        MaxValue = Math.Max(MaxValue, v)
                    Next

                End If

                'create three categories using 3 defined colors in options (color 0 is for undefined)

                Dim ShowLegendItems As Boolean = True

                Dim NumEntries As Integer = NumCats
                Dim debugtext As String = ""

                For i = 0 To NumEntries - 1

                    'create new category for symbology
                    Dim range As Double = (MaxValue - MinValue) / NumEntries
                    Dim start As Double = MinValue + range * i
                    Dim [end] As Double = start + range

                    Dim fillcolor As Color = GISForm.GradientColor(i * 1.0 / (NumEntries - 1))
                    Dim pc As New PolygonCategory(fillcolor, Color.Transparent, 0)
                    If i = NumEntries - 1 Then
                        pc.FilterExpression = $"[Result]>={start.ToString(Globalization.CultureInfo.InvariantCulture)}"
                    Else
                        pc.FilterExpression = $"[Result]>={start.ToString(Globalization.CultureInfo.InvariantCulture)} and [Result]<{[end].ToString(Globalization.CultureInfo.InvariantCulture)}"
                    End If

                    debugtext &= pc.FilterExpression & vbCrLf

                    pc.LegendText = String.Format("{0:N2} - {1:N2}", start, [end])
                    pc.LegendItemVisible = ShowLegendItems

                    'add category to scheme
                    ps.AddCategory(pc)

                Next

                Dim pcNull As New PolygonCategory(Color.Gray, Color.DarkGray, 1)
                pcNull.FilterExpression = $"[Result] IS NULL"
                pcNull.LegendText = "<Null>"
                pcNull.LegendItemVisible = ShowLegendItems
                ps.AddCategory(pcNull)

                ps.ResumeChangeEvent()
                ps.ResumeEvents()

                'set scheme to layer
                Try
                    lyrGrid.Symbology = ps
                Catch ex As Exception
                    lyrGrid.Symbology = Nothing
                    Throw New Exception("Unable to create layer scheme:" & vbCrLf & debugtext, ex)
                End Try

                GISForm.GISMap.Refresh()
                Debug.Print("CreateScheme: " & Now.Subtract(starttime).TotalSeconds)

            End With

        Catch ex As Exception
            ErrorMsg(, ex)
        End Try
    End Sub

    ''' <summary>
    ''' Uses the Douglas Peucker algorithm to reduce the number of points using DotSpatial data types. Note that interior holes are NOT preserved
    ''' </summary>
    ''' <param name="Poly">The polygon.</param>
    ''' <param name="Tolerance">The tolerance.</param>
    Private Shared Function DouglasPeuckerReduction(Poly As Polygon, Tolerance As Double) As Polygon
        Dim lstPoints1 As New List(Of PointF)
        For Each pt In Poly.ExteriorRing.Coordinates
            lstPoints1.Add(New PointF(pt.X, pt.Y))
        Next
        If lstPoints1(0).Equals(lstPoints1.Last) Then lstPoints1.RemoveAt(lstPoints1.Count - 1) 'remove last point if same as first (algorithm does not handle correctly)
        Dim lstPoints2 = DouglasPeuckerReduction(lstPoints1, Tolerance)
        Dim lstCoords2 As New List(Of Coordinate)
        For Each pt In lstPoints2
            lstCoords2.Add(New Coordinate(pt.X, pt.Y))
        Next
        Dim poly2 As New Polygon(lstCoords2) 'don't need to close up polygon
        Return poly2
    End Function

    ''' <summary>
    ''' Uses the Douglas Peucker algorithm to reduce the number of points.
    ''' Taken from http://www.codeproject.com/Articles/18936/A-Csharp-Implementation-of-Douglas-Peucker-Line-Ap
    ''' </summary>
    ''' <param name="Points">The points.</param>
    ''' <param name="Tolerance">The tolerance.</param>
    Private Shared Function DouglasPeuckerReduction(Points As List(Of PointF), Tolerance As Double) As List(Of PointF)
        If Points Is Nothing OrElse Points.Count < 3 Then
            Return Points
        End If

        Dim firstPoint As Int32 = 0
        Dim lastPoint As Int32 = Points.Count - 1
        Dim pointIndexsToKeep As New List(Of Int32)()

        'Add the first and last index to the keepers
        pointIndexsToKeep.Add(firstPoint)
        pointIndexsToKeep.Add(lastPoint)

        'The first and the last point cannot be the same
        While Points(firstPoint).Equals(Points(lastPoint))
            lastPoint -= 1
        End While

        DouglasPeuckerReduction(Points, firstPoint, lastPoint, Tolerance, pointIndexsToKeep)

        Dim returnPoints As New List(Of PointF)()
        pointIndexsToKeep.Sort()
        For Each index As Int32 In pointIndexsToKeep
            returnPoints.Add(Points(index))
        Next

        Return returnPoints
    End Function

    ''' <param name="points">The points.</param>
    ''' <param name="firstPoint">The first point.</param>
    ''' <param name="lastPoint">The last point.</param>
    ''' <param name="tolerance">The tolerance.</param>
    ''' <param name="pointIndexsToKeep">The point index to keep.</param>
    Private Shared Sub DouglasPeuckerReduction(points As List(Of PointF), firstPoint As Int32, lastPoint As Int32, tolerance As Double, ByRef pointIndexsToKeep As List(Of Int32))
        Dim maxDistance As [Double] = 0
        Dim indexFarthest As Int32 = 0

        For index As Int32 = firstPoint To lastPoint - 1
            Dim distance As [Double] = PerpendicularDistance(points(firstPoint), points(lastPoint), points(index))
            If distance > maxDistance Then
                maxDistance = distance
                indexFarthest = index
            End If
        Next

        If maxDistance > tolerance AndAlso indexFarthest <> 0 Then
            'Add the largest point that exceeds the tolerance
            pointIndexsToKeep.Add(indexFarthest)

            DouglasPeuckerReduction(points, firstPoint, indexFarthest, tolerance, pointIndexsToKeep)
            DouglasPeuckerReduction(points, indexFarthest, lastPoint, tolerance, pointIndexsToKeep)
        End If
    End Sub

    ''' 
    ''' The distance of a point from a line made from point1 and point2.
    ''' <param name="Point1">The PT1.</param>
    ''' <param name="Point2">The PT2.</param>
    ''' <param name="Point">The p.</param>
    ''' <returns></returns>
    Public Shared Function PerpendicularDistance(Point1 As PointF, Point2 As PointF, Point As PointF) As Double
        'Area = |(1/2)(x1y2 + x2y3 + x3y1 - x2y1 - x3y2 - x1y3)|   *Area of triangle
        'Base = v((x1-x2)²+(x1-x2)²)                               *Base of Triangle*
        'Area = .5*Base*H                                          *Solve for height
        'Height = Area/.5/Base

        Dim area As [Double] = Math.Abs(0.5 * (Point1.X * Point2.Y + Point2.X * Point.Y + Point.X * Point1.Y - Point2.X * Point1.Y - Point.X * Point2.Y - Point1.X * Point.Y))
        Dim bottom As [Double] = Math.Sqrt(Math.Pow(Point1.X - Point2.X, 2) + Math.Pow(Point1.Y - Point2.Y, 2))
        Dim height As [Double] = area / bottom * 2

        Return height

    End Function

End Class
